#ifndef _SYSTEM_GAMMOW_H_
#define _SYSTEM_GAMMOW_H_
#include<string>
#include<fstream>
#include"spsolver.h"
#include"system_ho.h"
class System_Tz_Gammow:public System_Tz<SP_K,KDataType>
{
 public:
  typedef SP_K OrbitalType;
  typedef System_Tz_HO System_TzType;
  const static int MatEledim;//< dimension of TwoBodyMatEle.
  typedef TwoBodyMatElement<KDataType> TwoBodyMatEle;
  typedef Matrix<TwoBodyMatEle,Dynamic,Dynamic> Mat2B;

 System_Tz_Gammow():Channels(spsolver.Channels){}  
 System_Tz_Gammow(const System_TzType * Ptr,int _A,int _Z):pSystem_Tz(Ptr),spsolver(Ptr,_A,_Z),Channels(spsolver.Channels),A(_A),Z(_Z),realA(_A),hbar_omega(Ptr->hbar_omega){}

  void setupOrbitals();
  void setupTwoBodyMat();
  void setAZ(int _A,int _Z,int _realA=-1)
  {
    if(_realA==-1) realA=_A;
    else realA=_realA;
    A=_A;
    Z=_Z;
    spsolver.setAZ(_A,_Z,realA);
  }

  void setup(const System_TzType * Ptr,int _A,int _Z)
  {
    pSystem_Tz=Ptr;
    hbar_omega=Ptr->hbar_omega;
    A=_A;
    Z=_Z;
    realA=_A;
    spsolver.setup(Ptr,_A,_Z);
    setupOrbitals();
    setupTwoBodyStates();
    setupTwoBodyMat();
  }
    
  void setup()
  {
    setupOrbitals();
    setupTwoBodyStates();
    //    setupTwoBodyMat();
  }
  ///interface for method beyond HF.
  KDataType get1B(int bra,int ket) const;
  KDataType get2B(int a,int b,int c,int d,int J) const
  {
    TwoBodyMatEle MatEle;
    //get2BmatAt(a,b,c,d,J,MatEle);
    get2BmatOnNeed(a,b,c,d,J,MatEle);
    if(a==b) MatEle*=sqrt(2.0);
    if(c==d) MatEle*=sqrt(2.0);
    return MatEle[0];
  }

  KDataType get1B_p2(int bra,int ket) const;

  void get2BmatOnNeed(int a,int b,int c,int d,int J,TwoBodyMatEle & MatEle) const;



  void printOrbitals(int num=-1) const;
  
  void print1B(const std::string outfile) const;
  void print2B(const std::string outfile) const;
  /* //maybe usefull operator */
  /* ///only radial part */
  /* KDataType rL(int bra,int ket,int L) const; */
  /* KDataType rms() const */
  /* { */
  /*   KDataType val=0; */
  /*   for(int i=0;i<=FermiSurface;i++) */
  /*     { */
  /* 	val+=rL(i,i,2); */
  /*     } */
  /*   return val/A; */
  /* } */
  /*   /\* */
  /*   Q0=sum_i{r_i^2},S(E0)=2\hbar^2/m*A*<r^2> */
  /* *\/ */
  /* ///be careful with the convention, Q0 is r^2, while QL proportion to r^L YL */
  /* KDataType Q(int p,int h,int J,int isospin) const; */
  /* KDataType Q1_v(int p,int h) const; */
  /* KDataType Q1_s(int p,int h,KDataType Rms) const; */


  double hbar_omega;
  
  int A,Z;
  int realA;
  int FermiSurface;
  KDataType HFenergy;
  const System_TzType * pSystem_Tz;
  SPSolver spsolver;
  vector<SPChannel> & Channels;
};


#endif
